1 



Generalized Rate-Code Model for Neuron Ensembles 
with Finite Populations 

Hideo Hasegawa0 

Department of Physics, Tokyo Gakugei University 
Koganei, Tokyo 184-8501, Japan 

(February 6, 2008) 

Abstract 

We have proposed a generalized Langevin-type rate-code model subjected 
to multiplicative noise, in order to study stationary and dynamical properties 
of an ensemble containing finite N neurons. Calculations using the Fokker- 
Planck equation (FPE) have shown that owing to the multiplicative noise, our 
rate model yields various kinds of stationary non-Gaussian distributions such 
as gamma, inverse-Gaussian-like and log-normal-like distributions, which have 
been experimentally observed. Dynamical properties of the rate model have 
been studied with the use of the augmented moment method ( AMM) , which 
was previously proposed by the author with a macroscopic point of view for 
finite-unit stochastic systems. In the AMM, original A^-dimensional stochastic 
differential equations (DEs) are transformed into three-dimensional determin- 
istic DEs for means and fluctuations of local and global variables. Dynamical 
responses of the neuron ensemble to pulse and sinusoidal inputs calculated by 
the AMM are in good agreement with those obtained by direct simulation. 
The synchronization in the neuronal ensemble is discussed. Variabilities of 
the firing rate and of the interspike interval (ISI) are shown to increase with 
increasing the magnitude of multiplicative noise, which may be a conceivable 
origin of the observed large variability in cortical neurons. 
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1 Introduction 



Neurons in a brain communicate information, emitting spikes which propagate 
through axons and dendrites to neurons at the next stage. It has been a long- 
standing controversy whether information of neurons is encoded in the firing rates 
(rate code) or in the more precise firing times (temporal code) [TJ, [2j [3] . Some experi- 
mental results having been reported seem to support the former code while some the 
latter [HO [3]. In particular, a recent success in brain-machine interface (BMI) [I] [5] 
suggests that the population rate code is employed in sensory and motor neurons 
while it is still not clear which code is adopted in higher-level cortical neurons. 

Experimental observations have shown that in many areas of the brain, neurons 
are organized into groups of cells such as columns in the visual cortex |6|. A small 
patch in cortex contains thousands of similar neurons, which receive inputs from the 
same patch and other patches. There are many theoretical studies on the property 
of neuronal ensembles consisting of equivalent neurons, with the use of spiking neu- 
ron models or rate-code models (for a review on neuronal models, see [7J; related 
references therein). In the spiking neuron model, the dynamics of the membrane 
potential of a neuron in the ensemble is described by the Hodgkin-Huxley (HH)- 
type nonlinear differential equations (DEs) [5] which express the conductance-based 
mechanism for firings. Reduced, simplified models such as the integrate-and-fire 
(IF) and FitzHugh-Nagumo (FN) models have been also employed. In contrast, in 
the rate-code model, neurons are regarded as transducers between input and output 
signals, both of which are expressed in terms of spiking rates. 

Computational neuroscientists have tried to understand the property of ensem- 
ble neurons by using the two approaches: direct simulations (DSs) and analytical 
approaches. DS calculations have been performed for large-scale networks mostly 
described by the simplest IF model. Since the computational time of DS grows as 
N 2 with N, the size of the ensemble, a large-scale DS with more realistic models be- 
comes difficult. Although DS calculations provide us with useful insight to the firing 
activity of the ensemble, it is desirable to have results obtained by using analytical 
approaches. 

Analytical or semi-analytical calculation methods for neuronal ensembles have 
been proposed by using the mean- field (MF) method 0, [10l [11], the population- 
density approaches p2]-[l6], the moment method [17] and the augmented moment 
method (AMM) [TS] (details of the AMM will be discussed shortly). It is interesting 
to analytically obtain the information of the firing rate or the interspike interval 
(ISI), starting from the spiking neuron model. It has been shown that the dynamics 
of the spiking neuron ensemble may be described by DEs of a macroscopic variable 
for the population density or spike activity, which determines the firing rate of 
ensemble neurons [12]- [IS]- By using the f-I relation between the applied dc current 
I and the frequency / of autonomous firings, the rate-code model for conduction- 
based neuron models is derived [19j [201 [21]. When we apply the Fokker-Planck 
equation (FPE) method to the neuron ensemble described by the IF model, the 
averaged firing rate R(t) is expressed by P(V, 9, t) which denotes the distribution 
probability of the averaged membrane potential V with the threshold 6 for the firing 
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[22J 



It is well known that neurons in brains are subjected to various kinds of noise, 
though their precise origins are not well understood. The response of neurons to 
stimuli is expected to be modified by noise in various ways. Indeed, although firings 
of a single in vitro neuron are reported to be precise and reliable [23], those of in 
vivo neurons are quite unreliable due to noisy environment. The strong criticism 
against the temporal code is that it is not robust against noise, while the rate code 
is robust. 

It is commonly assumed that there are two types of noise: additive and multi- 
plicative noise. The magnitude of the former is independent of the state of variable 
while that of the latter depends on its state. Interesting phenomena caused by 
the two noise have been investigated [24j. It has been found that the property of 
multiplicative noise is different from that of additive noise in some respects. (1) Mul- 
tiplicative noise induces a phase transition, creating an ordered state, while additive 
noise works to destroy the ordering [25] [2E]- (2) Although the probability distribu- 
tion in stochastic systems subjected to additive white noise follows a Gaussian, mul- 
tiplicative white noise generally yields a non-Gaussian distribution [27]-[32]. (3) The 
scaling relation of the effective strength for additive noise given by (3(N) = (3(1)/ y/~N 
is not applicable to that for multiplicative noise: a(N) ^ a(l)/\/~N, where a(N) 
and /3(N) denote effective strengths of multiplicative and additive noise, respec- 
tively, in the iV-unit system [33]. A naive approximation of the scaling relation for 
multiplicative noise: a(N) = ac(l)/y/N as adopted in Ref. [22], yields the result 
which does not agree with that of DS [33J. 

In this paper, we will study the property of neuronal ensembles based on the 
rate-code hypothesis. Rate models having been proposed so far, are mainly given 



where Ti(t) (> 0) denotes a firing rate of a neuron i [i = 1 to N), A the relaxation 
rate, Wij the coupling strength, H(x) the gain function, Jj(t) an external input, 
and (3 expresses the magnitude of additive white noise of ^(t) with the correlation: 
< £■,-(£') >= 5ij5(t — t'). The rate model as given by Eq. (1) has been adopted 
in many models based on neuronal population dynamics. The typical rate model is 
the Wilson-Cowan model, with which the stability of an ensemble consisting of exci- 
tatory and inhibitory neurons is investigated [12j[34|. The rate model given by Eq. 
(1) with H(x) = x is the Hopfield model [35J, which has been extensively adopted 
for studies on the memory in the brain incorporating the plasticity of synapses into 
Wij. DS calculations have been performed, for example, for a study of the population 
coding for iV = 100 |4J. Analytical studies of Eq. (1) are conventionally made for the 
case of N = oo, adopting the FPE method with MF and diffusion approximations. 
The stationary distribution obtained by the FPE for Eq. (1) generally follows the 
Gaussian distribution. 

ISI data obtained from experiments have been fitted by a superposition of some 
known probability densities such as the gamma, inverse- Gaussian and log-normal 
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distributions [36]-[l0]. The gamma distribution with parameters A and /i is given 
by 

W^H^a^exp^-^), (2) 

which is derived from a simple stochastic IF model with additive noise for Poisson 
inputs [37], T(x) being the gamma function. For A = 1 in Eq. (2), we get the ex- 
ponential distribution describing a standard Poisson process. The inverse Gaussian 
distribution with parameters A and /i given by 



Pig(x) 




x(x — ny 



(3) 



is obtained from a stochastic IF model in which the membrane potential is repre- 
sented as a random walk with drift [36J . The log- normal distribution with parameters 
fj, and a given by 



Pln{x) = rj— exp 



V2 



7ia 2 x 



(logx — 
2^2 



(4) 



is adopted when the log of ISI is assumed to follow a Gaussian form [38] . Fittings 
of experimental ISI data to a superposition of these probability densities have been 
extensively discussed in the literature [36]-[4"0]. 

The purpose of the present paper is to propose and study the generalized, phe- 
nomenological rate model [Eqs. (5) and (6)]. We will discuss ensembles with finite 
populations, contrary to most existing analytical theories except some ones {e.g. 
Ref. [15]), which discuss ensembles with infinite N . The stationary distribution of 
our rate model will be discussed by using the FPE method. It is shown that owing to 
the introduced multiplicative noise, our rate model yields not only the Gaussian dis- 
tribution but also non-Gaussian distributions such as gamma, inverse-Gaussian-like 
and log-normal-like distributions. 

The dynamical properties of our rate model will be studied by using the aug- 
mented moment method (AMM) which was previously proposed by the present 
author [TBI 1331 HJJ. Based on a macroscopic point of view, Hasegawa [18] has pro- 
posed the AMM, which emphasizes not the property of individual neurons but rather 
that of ensemble neurons. In the AMM, the state of finite iV-unit stochastic ensem- 
bles is described by a fairly small number of variables: averages and fluctuations 
of local and global variables. For example, the number of deterministic equation in 
the AMM becomes three for a iV-unit Langevin model. The AMM has been suc- 
cessfully applied to a study on the dynamics of the Langevin model and stochastic 
spiking neuron models such as FN and HH models, with global, local or small-world 
couplings (with and without transmission delays) [4"2"]-[4"6]. 

The AMM in [18] was originally developed by expanding variables around their 
stable mean values in order to obtain the second-order moments both for local and 
global variables in stochastic systems. In recent papers [331 E] > we have reformu- 
lated the AMM with the use of the FPE to discuss stochastic systems subjected to 
multiplicative noise: the FPE is adopted to avoid the difficulty due to the Ito versus 
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Stratonovich calculus inherent to multiplicative noise. In the present paper, a study 
on the Langevin model with multiplicative noise made in [41], has been applied to 
an investigation on the firing properties of neuronal ensembles. Our method aims 
at the same purpose to effectively study the property of neuronal ensembles as the 
approaches developed in Refs. [T2"]-[T6][r9]-[2~Tj. 

The paper is organized as follows. In Sec. 2, we discuss the generalized rate model 
for an ensemble containing N neurons, investigating its stationary and dynamical 
properties. Some discussions are presented in Sec. 3, where variabilities of the rate 
and ISI are calculated. The final Sec. 4 is devoted to our conclusion. 



2 Property of neuron ensembles 
2.1 Generalized rate-code model 

For a study of the property of a neuron ensemble containing finite iV neurons, we 
have assumed that the dynamics of the firing rate r«(t) (> 0) of a neuron i (i — 1 
to N) is given by 



dt 
with 



(IT • 

1 -- F( ri ) + H( Ul )+aG(r t ) m (t)+f3at), (5) 



«*(*) = f|) E rM+Ut). (6) 



Here F(x), G(x) and H(x) are arbitrary functions of x, Z (= N — 1) denotes the co- 
ordination number, I^(t) an input from external sources and w the coupling strength: 
a and f3 express the strengths of additive and multiplicative noise, respectively, given 
by and i]i(t) expressing zero- mean Gaussian white noise with correlations given 
by 

< m(t) Vj(t') > = Sydit-t), (7) 
<&(*)6(0> = (8) 

<ifc(t)fc(f)> = o. (9) 

The rate model in Eq. (1) adopts F(x) = —Xx and G(x) = (no multiplicative 
noise) . 

The gain function H (x) expresses the response of the firing rate (r») to a synaptic 
input field (uj). It has been theoretically shown in [47] that when spike inputs with 
the mean ISI (Tj n ) are applied to an HH neuron, the mean ISI of output signals 
(T out ) is T out = T in for T in ~ 15 ms and T out ~ 15 ms for T in ~ 15 ms. This is 
consistent with the recent calculation for HH neuron multilayers |48j , which shows a 
nearly linear relationship between the input (r» n ) and output rates (r out ) for r in < 60 
Hz (Fig. 3 of Ref. [H]). It is interesting that the r in - r out relation is continuous 
despite the fact that the f-I relation of the HH neuron shows a discontinuous, type- 
II behavior according to Ref. [8]. In the literature, two types of expressions for 
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H(x) have been adopted so far. In the first category, sigmoid functions such as 
H(x) = 1/(1 + e~ x ) (e.g. [12]) and arctan(x) (e.g. [19]) have been adopted. In 
the second category, gain functions such as H(x) = (x — x c )Q(x — x c ) (e.g. |21j) 
have been employed, modeling the f-I function for the frequency / of autonomous 
oscillation against the applied dc current J, x c expressing the critical value and 
<d(x) the Heaviside function: Q(x) = 1 for x > and otherwise. The nonlinear, 
saturating behavior in H(x) arises from the property of the refractory period (r r ) 
because spike outputs are prevented for tf < t < tf + r r after firing at t — tf. In 
this paper, we have adopted a simple expression given by [50J 



H(x) 



x 



(10) 



Vx^TT' 

although our results to be presented in the following sections, are expected to be 
valid for any choice of H(x). 



2.2 Stationary properties 
2.2.1 Distribution of r 

By employing the FPE, we may discuss the stationary distribution p(r) for w = 
and Ii(i) = I, which is given by [30], [31] 

~a 2 G(r) 2 P 2 ' 



\np(r) oc X(r) + Y(r) - I 1 - ^ j In 



with 



X(r) 
Y(r) 



2 j dr 
2 I dr 



F(r) 



a 2 G(r) 2 + p 2 

H(D 

a 2 G(r) 2 + p 2 
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12) 



(13) 



where = and 1 for Ito and Stratonovich representations, respectively. Hereafter 
we mainly adopt the Stratonovich representation. 



Case I F(x) = —Xx and G(x) = x 

For the linear Langevin model, we get 



p(r) oc 



1 + 



' a 2 r 2S 



2 



with 



-(A/a 2 +l/2) 



Y(r) = (¥L\ arctan 



(14) 



(15) 



where H = H(I). In the case of H = Y(r) = 0, we get the g-Gaussian given by 

EOT 



p(r) oc 1 — (1 — q)^r* 



l-q 



(16) 
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with 



7 
<1 

We examine some limiting cases of Eq. (14) as follows. 

(a) For a = and (3 ^ (i.e. additive noise only), Eq. (14) yields 

p(r) oc e -^ (r ~ H/X)2 . (19) 

(b) For (3 — and a^O (i.e. multiplicative noise only), Eq. (14) leads to 

P (r) OC r -(2A/« 2 +D e -(2^)/r_ (2Q) 

Distributions p(r) calculated with the use of Eqs. (14)-(20) are plotted in Figs. 
l(a)-l(c). The distribution p(r) for a = 0.0 (without multiplicative noise) in Fig. 
1(a) shows the Gaussian distribution which is shifted by an applied input / = 0.1. 
When multiplicative noise is added (a ^ 0), the form of p(r) is changed to the 
g-Gaussian given by Eq. (16). Fig. 1(b) shows that when the magnitude of additive 
noise (3 is increased, the width of p(r) is increased. Fig. 1(c) shows that when the 
magnitude of external input I is increased, p(r) is much shifted and widely spread. 
Note that for a = 0.0 and (3^0 (additive noise only), p(r) is simply shifted without 
a change in its shape when increasing / [Eq. (19)]. 



2A + Q 2 
2(3 2 ' 
2A + 3a 2 
2A + o< 2 ' 



(17) 



Case II F(x) = -\x a and G(x) = x b (a, b > 0) 

The special case of a = 1 and 6 = 1 has been discussed in the preceding case I 
[Eqs. (14)-(20)]. For arbitrary a (> 0) and b (> 0), the probability distribution p(r) 
given by Eqs. (11)-(13) becomes 



p(r) oc 



1 + 



:2b 



-1/2 



,X(r)+Y(r) 



with 



X(r) 
Y(r) 



/ 2Ar a+1 
~{(3 2 (a + l) 
(2Hr\ 



( a + 1 a + 1 



a 2 r 2bS 



2b ' 2b 



F V>&k +1 > 



a 2 r 2bS 



i3 2 



(21) 



(22) 



(23) 



where F(a, b, c; z) is the hypergeometric function. Some limiting cases of Eqs. (21)- 
(23) are shown in the following. 

(a) The case of H = Y(r) = was previously studied in 

(b) For a = and (3 ^ (i.e. additive noise only), we get 



p(r) oc exp 



2A 



f 2 (a + l) 



r a+1 + 



'2lf 



(24) 
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(c) For (3 = and a^O (i.e. multiplicative noise only), we get 

/ 2A \ / 2H 



p(r) oc r exp 



„a-2ft+l 



,-26+1 



oc r-( 2A / Q2 + fe )exp 



a 2 (a-26+l)/ V" 2 ^- 1 ) 
for a - 26 + 1 ^ 0, 26 - 1 ^ 
/ 2# 



-26+1 



a 2 (2b- 1 
/ 2A \ 



(25) 

, for a - 2b + 1 = (26) 



a 2 a 



for 26 - 1 = 



oc r 



-[2(A-//)/a 2 +l/2] 



for a - 26 + 1 = 0, 26 - 1 = 



(d) In the case of a = 1 and 6 = 1/2, we get 



/ (2A/3 2 /a 4 +2ff/a 2 -l/2) 

p(r) oc ( r + ^J 



cxp 



' 2A v 



cr 



which reduces, in the limit of a = 0, to 



p(r) oc exp 



H 



for a = 



(27) 
(28) 

(29) 
(30) 



Case III F(x) = — Xlnx and G(x) = x 1 / 2 
We get 



p(r) oc r exp 









( lnr -f) 
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for /3 = 



(31) 



Figure 2(a) shows distributions p(r) for case II and various a with fixed values 
of A = 1.0, 6 = 1.0, / = 0.1, a = 1.0 and (3 = 0.0 (multiplicative noise only). With 
more decreasing a, a peak of p(r) at r ~ 0.1 becomes sharper. Fig. 2(c) shows 
distributions p(r) for case II and various 6 with fixed values of A = 1.0, a = 1.0, 
/ = 0.1, a = 1.0 and (3 = 0.0 (multiplicative noise only). We note that a change in 
the 6 value yields considerable changes in shapes of p(r). Figs. 2(b) and 2(d) will 
be discussed shortly. 



2.2.2 Distribution of T 

When the temporal ISI T is simply defined by T = 1/r, its distribution n(T) is 
given by 



ir(T)=p 



1\ 1 



T J T 2 ' 



(32) 



We get various distributions of vr(T) depending on functional forms of F(x) and 
G(x). For F(x) = -\x, G(x) = x and (3 = 0, Eq. (26) yields 



tt(T) oc T^/" 2 - 1 ) exp 
8 



(2H 



T 



(33) 



which expresses the gamma distribution [Eq. (2)] [37J [29]. For F(x) = —Ax 2 , 
G(x) = x and (3 = 0, Eq. (25) leads to 



which is similar to the inverse Gaussian distribution [Eq. (3)] [35] • For F(x) = 
-Xlnx, G(x) = x 1 ' 2 and (3 = 0, Eq. (31) yields 



which is similar to the log-normal distribution [Eq. (4)] [38J. 

Figs. 2(b) and 2(d) show 7r(T) obtained from p(r) shown in Figs. 2(a) and 2(c), 
respectively, by a change of variable with Eq. (32). Fig. 2(b) shows that with more 
increasing a, the peak of tt(T) becomes sharper and moves left. We note in Fig. 
2(d) that the form of vr(T) significantly varied by changing b in G(x) = x b . 

2.2.3 Distribution of R 

When we consider the global variable R(t) defined by 




(34) 




(35) 



R(t) = ^ !>«(*), 



(36) 



the distribution P(R, t) for R is given by 




(37) 



Analytic expressions of P(R) are obtainable only for limited cases, 
(a) For (3 ^ and a = 0, P(R) is given by 




(38) 



where H = H(I). 

(b) For H = 0, we get @T] 




(39) 



with 




(40) 



where (f>(k) is the characteristic function for p(r) given by [ST] 




(41) 



2 



(A- | k |) 

r(iz) 



K v {\' | k |) 



(42) 
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with 



v = 4 ( 43 ) 

A' = K (44) 
a 

K v {x) expressing the modified Bessel function. 

Some numerical examples of P(R) are plotted in Figs. 3, 4 and 5. Figures 3(a) 
and 3(b) show P(R) for a = 0.0 and a = 0.5, respectively, when N is changed. For 
a = 0.0, P(R) is the Gaussian distribution whose width is narrowed by a factor of 
1/yN with increasing N. In contrast, P(R) for a = 0.5 is non-Gaussian, whose 
shape seems to approach a Gaussian for increasing N. These are consistent with 
the central-limit theorem. 

Effects of an external input I on p(r) and P(R) are examined in Figs. 4(a) and 
4(b). Figure 4(a) shows that in the case of a = 0.0 (additive noise only), p(r) and 
P(R) are simply shifted by a change in J. This is not the case for a ^ 0.0, for which 
p{r) and P(R) are shifted and widened with increasing J, as shown in Fig. 4(b). 

Figures 5(a) and 5(b) show effects of the coupling w on p(r) and P(R)- For 
a = 0.0, p{r) and P(R) are changed only slightly with increasing w. On the contrary, 
for a = 0.5, an introduction of the coupling significantly modifies p(r) and P(R) as 
shown in Fig. 5(b). 

2.3 Dynamical properties 

2.3.1 Augmented Moment Method (AMM) 

Next we will discuss the dynamical properties of the rate model by using the AMM 
[T8"t I3"3"t HTJ. By employing the FPE, we obtain equations of motion for moments: 
(rj), (rj rj), and (R 2 ) where R = (1/N) J2i r i- Then we get equations of motion for 
the three quantities of /i, 7 and p defined by [HI [33], [JT] 

V- = («> = ^E<n>, (45) 

% 

7 = ^E((^-/^) 2 ) ; (46) 

i 

P = ((R-P) 2 ), (47) 

where \x expresses the mean, 7 the averaged fluctuations in local variables (r^) and 
p fluctuations in the global variable (R). We get (for details see [33[ |4T] ) 

^ = fo + hi + h + [9o9i + 3(^i^2 + 9093)7], (48) 

d 7 nf , oh fwN\ ( 7 

+ (<P+l)(g 2 1 + 2g g 2 )a 2 1 + a 2 g 2 +(3 2 , (49) 
^ = 2f l p + 2h l wp+(<j ) +l)(g 2 + 2g g 2 )a 2 p + -^(a 2 g 2 + p 2 ) J (50) 
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where f t = (l/£\){d e F{p)/dx e ), g e = {l/£\){d e G{p)/dx e ), h e = (l/^!)(^(u)/9u<) 
and u = wpi + I. Original iV-dimensional stochastic DEs given by Eqs. (5) and (6) 
are transformed to the three-dimensional deterministic DEs given by Eqs. (48)- (50). 
When we adopt 

F(x) = -Ax, (51) 
G{x) = x, (52) 

Eqs. (48)- (50) are expressed in the Stratonovich representation (0 = 1) by 

- = -X li + ho + —, (53) 

% = -2A 7 + ^(p-^) + 2^ 7 + aV + ^ (54) 

f f = -2\p + 2h 1 wp + 2a 2 p + ^f + ^, (55) 

where h = u/^u 2 + 1, h x = l/(u 2 + if' 2 and h 2 = -(3m/2)/(m 2 + if' 2 with 
u = wp + /. 

Before discussing the dynamical properties, we study the stationary properties 
of Eqs. (53)-(55). We get the stationary solution given by 

" = (A-aV2)' (56) 



7 = 



2(A - a 2 + whjZ) 

' (aV + /3 2 ) 
2iV(A -a 2 - whxY 



1 + ^ 



Z(A — a 2 — u»/ii) 



(57) 
(58) 



where Eq. (56) expresses the fifth-order algebraic equation of /i. The stability of Eqs. 
(53)-(55) around the stationary solution may be shown by calculating eigenvalues 
of their Jacobian matrix, although actual calculations are tedious. 

Figure 6 shows the N dependences of 7 and p in the stationary state for four 
sets of parameters: (a,(3,w) = (0.0,0.1,0.0) (solid curves), (0.5, 0.1, 0.0) (dashed 
curves), (0.0, 0.1, 0.5) (chain curves) and (0.5, 0.1, 0.5) (double-chain curves), with 
(3 = 0.1, A = 1.0 and I = 0.1. For all the cases, p is proportional to iV -1 , which is 
easily seen in Eq. (58). In contrast, 7 shows a weak iV dependence for a small iV 
(< 10). It is noted that y/j and yfp approximately express the widths of p(r) and 
P(R), respectively. The iV-dependence of p in Fig. 6 is consistent with the result 
shown in Figs. 3(a) and 3(b), and with the central-limit theorem. 



2.3.2 Response to pulse inputs 

We have studied the dynamical properties of the rate model, by applying a pulse 
input of I = I(t) given by 

i(t) = Ae(t-t 1 )e(t 2 -t) + i<- b \ (59) 
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with A = 0.5, h = 40, t 2 = 50 and = 0.1 which expresses the background input. 
Figs. 7(a), 7(b) and 7(c) show the time dependences of p, 7 and p for F(x) = —\x 
and G(x) = x when the input pulse I(t) given by Eq. (59) is applied [52]: solid 
and dashed curves show the results of the AMM and DS averaged over 1000 trials, 
respectively, with a = 0.5, (3 = 1.0, w = 0.5 and N = 10. Figs. 7(b) and 7(c) show 
that an applied input pulse induces changes in 7 and p. This may be understood 
from 2a 2 terms in Eqs. (54) and (55). The results of the AMM shown by solid 
curves in Figs. 7(a)-(c) are in good agreement with DS results shown by dashed 
curves. Figure 7(d) will be discussed shortly. 

It is possible to discuss the synchrony in a neuronal ensemble with the use of 
7 and p defined by Eqs. (46) and (47) [15] . In order to quantitatively discuss the 
synchronization, we first consider the quantity given by 



When all neurons are firing with the same rate (the completely synchronous state), 
we get Tiit) = R(t) for all i, and then P(t) = in Eq. (60). On the contrary, we get 
P(t) = 2(1 - I/AO7 = P (t) in the asynchronous state where p = 7/JV [i~8ll4Tj . We 
may define the synchronization ratio given by [18] 



which is and 1 for completely asynchronous (P = Pq) and synchronous states 
(P = 0), respectively. Figure 7(d) shows the synchronization ratio S(t) for / ~f(t) and 
p(t) plotted in Figs. 7(b) and 7(c), respectively, with a = 0.5, (3 = 1.0, w = 0.5 and 
iV = 10. The synchronization ratio at t < 40 and t > 60 is 0.15, but it is decreased 
to 0.03 at 40 < t < 50 by an applied pulse. This is because by an applied pulse, 7 
is more increased than p and the ratio of p/7 is reduced. The synchronization ratio 
vanishes for w = 0, and it is increased with increasing the coupling strength [T8| |4"T]. 

Next we show some results for different indices of a and b in F(x) = —\x a and 
G(x) = x b . Fig. 8(a) shows the time dependence of p for (a, b) = (1, 1) (solid curve) 
and (a, b) = (2, 1) (dashed curve) with a = 0.0, = 0.1, w = 0.0 and N = 10. The 
saturated magnitude of p for a = 0.5 is larger than that for a = 0.0. Solid and 
dashed curves in Fig. 8(b) show p for (a, b) = (1,1) and (1,0.5), respectively, with 
a = 0.5, (3 = 0.001, iV = 10 and w = 0.0. Both results show similar responses to an 
applied pulse although p for a background input of 1^ = 0.1 for (a, b) = (1,0.5) is 
a little larger than that for (a, b) = (1, 1). 

2.3.3 Response to sinusoidal inputs 

We have applied also a sinusoidal input given by 



Pit) 



1 

iV 2 



£<[r*(t)-r,(t)] 2 >=2[ 7 (t)-p(f)]. 



(60) 




(61) 




(62) 
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for F(x) = -A x and G(x) = x with A = 1.0, A = 0.5, I® = 0.1, and T p = 10 and 
20. Time dependences of \x for T p = 20 and T p = 10 are plotted in Figs. 9(a) and 
9(b), respectively, with a = 0.5, (3 = 1.0, w = 0.0 and N = 10. AMM results of //(t) 
shown by solid curves in Figs. 9(a) and (b) are indistinguishable from DS results 
(with 100 trials) shown by dashed curves [52], chain curves denoting sinusoidal input 
I{t). As the period of T p becomes shorter, the magnitude of [i becomes smaller. The 
delay time of /i(t) against an input I(t) is about r d ~ 1.0 (= 1/A) for both T p = 10 
and T p = 20. 



3 Discussion 

We may calculate variabilities of r and T, by using their distributions of p(r) and 
n(T), which have been obtained in Sec. 2. For example, in the case of F(x) = —Xx 
and G(x) = x, the distribution of p(r) for (3 = 0.0, w = 0.0 and H = I given by Eq. 
(26) leads to 

< r > - (X^J2 Y (63) 
^ - « r -W = 2iX-JimX-a>y (64) 



c v = - = (65) 
( r ) ^2{X - a 2 ) 

The relevant gamma distribution for ISI, n(T), given by Eq. (33) yields 

(T) = j, (66) 

Xa 2 

(5T 2 ) = ((T-(T)) 2 ) = W , (67) 

(68) 



CO V2X 

Equations (65) and (68) show that both c v and c v are increased with increasing the 
magnitude (a) of multiplicative noise. 

It has been reported that spike train variability seems to correlated with loca- 
tion in the processing hierarchy [53] . A large value of c v is observed in hippocampus 
(c v ~ 3) [54j whereas c v is small in cortical neurons {c v ~ 0.5 — 1.0) and motor neu- 
rons (c v ~ 0.1) [55j[56j. In order to explain the observed large c v , several hypotheses 
have been proposed: (1) a balance between excitatory and inhibitory inputs [57j[58j, 
(2) correlated fluctuations in recurrent networks [59], (3) the active dendrite con- 
ductance [60], and (4) a slowly decreasing tail of input ISI of T~ d (d > 0) at large 
T [61] . Our calculation shows that multiplicative noise may be an alternative origin 
(or one of origins) of the observed large variability. We note that the variability 
of r is given by c v = in the AMM (e.g. Eq. (65) agrees with for /i 
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and 7 given by Eqs. (56) and (57), respectively, with w = (3 = 0). It would be 
interesting to make a more detailed study of the variability for general F(x) and 
G(x) as discussed in Sec. 2. 

We have proposed the generalized rate-code model given by Eqs. (5) and (6), in 
which the relaxation process is given by a single F(x). Instead, when the relaxation 
process consists of two terms: 

F(x)^c l F 1 (x)+c 2 F 2 (x), (69) 

with c\ + c 2 = 1, the distribution becomes 

p(r) = Mr)] Cl b 2 (r)] C2 , (70) 

where Pfc(r) (k = 1,2) denotes the distribution only with F(x) = Fi(x) or F(x) = 
Fi(x). In contrast, when multiplicative noise arises from two independent origins: 

axr](t) — > ciaixr]i(t) + c 2 a 2 xr] 2 (t), (71) 

the distribution for f3 = H = becomes 

p(r) oc r ~[ 2A /(cia?+ c 2 a i)+ 1 ], (72) 

Similarly, when additive noise arises from two independent origins: 

m)^aPMt) + c 2 (3 2 Ut), (73) 

the distribution for a = H = becomes 

p(r) a e - A/(cia ? +C2a ^. (74) 

Equations (70), (72) and (74) are quite different from the form given by 

p(r) = cipi(r) + c 2 p 2 (r), (75) 

which has been conventionally adopted for a fitting of the theoretical distribution 
to that obtained by experiments [36]-[40j. 

4 Conclusion 

We have proposed the generalized rate-code model [Eqs. (5) and (6)], whose proper- 
ties have been discussed by using the FPE and the AMM. The proposed rate model 
is a phenomeno logical one and has no biological basis. As discussed in Sec. 1, the 
conventional rate model given by Eq. (1) may be obtainable from a spiking neu- 
ron model when we adopt appropriate approximations to DEs derived by various 
approaches such as the population-density method p2]-[T6] and others [19]- [21]. It 
would be interesting to derive our rate model given by Eqs. (5) and (6), starting 
from a spiking neuron model. The proposed generalized rate model is useful in dis- 
cussing stationary and dynamical properties of neuronal ensembles. Indeed, our rate 
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model has an interesting property, yielding various types of stationary non-Gaussian 
distributions such as gamma, inverse- Gaussian and log-normal distributions, which 
have been experimentally observed [36]- [40]. It is well known that the Langevin- 
type model given by Eq. (1) cannot properly describe fast neuronal dynamics at the 
characteristic times r c shorter than r (= 1/A). This is, however, not a fatal defect 
because we may evaded it, by adopting an appropriate r value of r < r c /10 for a 
given neuronal ensemble with r c . Actually, the dynamical properties of an ensemble 
consisting of excitatory and inhibitory neurons has been successfully discussed with 
the use of the Langevin-type Wilson-Cowan model [12] [34] (for recent papers using 
the Wilson-Cowan model, see [62]: related references therein). One of the disad- 
vantages of the AMM is that its applicability is limited to the case of weak noise 
because it neglects contributions from higher moments. 
On the contrary, the AMM has following advantages: 

(i) the dynamical properties of an iV-unit neuronal ensemble may be easily studied 
by solving three-dimensional ordinary DEs [Eqs. (48)-(50)], in which three quantities 
of /i, 7 and p have clear physical meanings, 

(ii) analytic expressions for DEs provide us with physical insight without numerical 
calculations (e.g. the N dependence of p follows the central-limit theorem [Eq. (58)], 
and 

(hi) the synchronization of the ensemble may be discussed [Eq. (61)]. 
As for the item (i), note that we have to solve the iV-dimensional stochastic Langevin 
equations in DS, and the (2N + l)-dimensional partial DEs in the FPE. Then the 
AMM calculation is very much faster than DS: for example, for the calculation shown 
in Fig. 9(a), the ratio of the computation time of the AMM to that of DS becomes 
t amm /ins ~l/30 000 [63]. We hope that the proposed rate model may be adopted 
for a wide class of study on neuronal ensembles described by the Wilson- Cowan- type 
model [64] . 
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Figure 1: (a) Distributions p(r) of the (local) firing rate r for various a with A = 1.0, 
= 0.1, / = 0.1 and w = 0.0, (b) p{r) for various (3 with A = 1.0, a = 1.0, I = 0.1 
and w = 0.0, and (c) p(r) for various / with A = 1.0, a = 0.5, (3 = 0.1 and w = 0.0. 

Figure 2: (a) Distributions p{r) of the (local) firing rate r and (b) vr(T) of the 
ISI T for a = 0.8 (chain curves), a = 1.0 (solid curves), a = 1.5 (dotted curves) 
and a = 2.0 (dashed curves) with A = 1.0, b = 1.0, / = 0.1, a = 1.0 and (3 = 0.0 
(multiplicative noise only), (c) p(r) of the (local) firing rate r and (d) 7r(T) of the 
ISI T for b = 0.5 (dashed curves), b = 1.0 (solid curves), b = 1.5 (dotted curves) 
and b = 2.0 (chain curves) with A = 1.0, a = 1.0, / = 0.1, a = 1.0 and f3 = 0.0 
(multiplicative noise only): results for b = 1.5 and b = 2 should be multiplied by 
factors of 2 and 5, respectively. 

Figure 3: Distributions P{R) of the (global) firing rate R for (a) a = 0.0 and (b) 
a = 0.5, with N = 1, 10 and 100: A = 1.0, /3 = 0.1, w = 0.0 and I = 0.1. 

Figure 4: Distributions p(r) (dashed curves) and P{R) (solid curves) for (a) a = 0.0 
and (b) a = 0.5 with / = 0.1 and / = 0.2: iV = 10, A = 1.0, (5 = 0.1 and w = 0.0. 

Figure 5: Distributions p(r) (dashed curves) and P{R) (solid curves) for (a) a = 0.0 
and (b) a = 0.5 with w = 0.0 and w = 0.5: iV = 10, A = 1.0, (5 = 0.1 and / = 0.1. 

Figure 6: The N dependence of 7 and p in the stationary states for four sets of 
parameters: (a,(3,w) = (0.0,0.1,0.0) (solid curves), (0.5,0.1,0.0) (dashed curves), 
(0.0,0.1,0.5) (chain curves) and (0.5,0.1,0.5) (double-chain curves): A = 1.0, iV = 
10 and / = 0.1. 

Figure 7: (Color online) Time courses of (a) /i(t), (b) 7(t), (c) p(t) and (d) S(t) 
for a pulse input I(t) given by Eq. (59) with A = 1.0, a = 0.5, (3 = 0.1, iV = 10 
and w = 0.5, solid and chain curves denoting results of AMM and dashed curves 
expressing those of DS result with 1000 trials. 

Figure 8: (a) Response of ix{t) to input pulse I(t) given by Eq. (59) for (a, b) = (1, 1) 
(solid curve) and (a, b) = (2, 1) (dashed curve) with a = 0.0, (3 = 0.1, iV = 10 and 
A = 1.0. (b) Response of /i(t) to input pulse I(t) for (a, b) = (1,1) (solid curve) 
and (a, b) = (1,0.5) (dashed curve) with a = 0.5, (3 = 0.001, iV = 10, A = 1.0 and 
w = 0.0. 

Figure 9: (Color online) Response of p,(t) to sinusoidal input I(t) (chain curves) 
given by Eq. (62) for (a) T p = 20 and (b) T p = 10 with A = 0.5, A = 1.0, a = 0.5, 
(3 — 0.1, w — 0.0 and TV = 10 (a — 1 and 6=1), solid and chain curves denoting 
fj,(t) of AMM and dashed curves expressing those of DS result with 100 trials. 
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